Digital Quantum Simulation of the Holstein Model in Trapped Ions 
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We propose the implementation of the Holstein model by means of digital methods in a linear 
chain of trapped ions. We show how the simulation fidelity scales with the generation of phononic 
excitations. We propose a decomposition and a stepwise trapped-ion implementation of the Holstein 
Hamiltonian. Via numerical simulations, we study how the protocol is affected by realistic gates. 
Finally, we show how measurements of the size of the simulated polaron can be performed. 
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Quantum simulators [TJ [2] are promising tools for the 
deep comprehension of complex quantum dynamics. In 
a quantum simulator, the higher control on the simu- 
lating system can allow to reproduce and recover non- 
trivial quantum behaviors. Recently, a significant boost 
to the field of quantum simulations has been provided 
by the use of digital approximations in trapped-ion se- 
tups [3l H] , based on stroboscopic decompositions of uni- 
tary operators [HIS]. However, the digital simulation of 
coupled bosonic-fermionic systems, naturally described 
by unbounded Hamiltonians, has not been considered. 

Strongly correlated quantum many body systems rep- 
resent a challenge to both computational and analytic 
methods. Among them, correlated fermionic-bosonic 
models are of critical relevance. The importance of cor- 
relation between electrons and ion vibrations has been 
proven for a large number of condensed-matter sys- 
tems [TJ. Their role in high-temperature superconduc- 
tors, as fullerides and cuprates, is still debated [8HI0]. In 
solid state systems, the correlation between the presence 
of electrons in a lattice and deformations of the latter 
can result in the formation of polarons: electrons and 
phonons can no longer be considered as stand-alone par- 
ticles. Depending on the strength of the electron-phonon 
couplings, the cloud of lattice displacements surrounding 
the electron can have different sizes. For strong cou- 
plings, the electrons can be trapped, with remarkable 
changes of global properties [TJ . The Holstein model [T2] 
has been proved to naturally describe the strong cou- 
pling case. This model has been recently addressed 
by heavy numerical simulations |13j and classical ana- 
log simulations for a reduced number of sites [13] , Per- 
turbation methods based on the Lang-Firsov approxima- 
tions [15] . valid in the strong coupling limit, are known 
since long times. The dimensionality of the underlying 
lattice also raises critical features [TB]. While involving 
a lot of efforts, the full and complete comprehension of 
the electron-phonon correlations is still an open problem. 

Trapped-ion systems are among the most controllable 
quantum systems. They offer remarkable computational 
power to perform quantum simulations exponentially 
faster than their classical counterparts [171132] . 



In this Letter, we propose the implementation of the 
Holstein Hamiltonian in a chain of trapped ions, using 
digital approximation methods. First, we address the 
problem of simulating with digital protocols unbounded 
Hamiltonians. Then, we provide a convenient decomposi- 
tion of the Holstein Hamiltonian, in that each step can be 
implemented in a trapped-ion setup. We discuss a possi- 
ble experimental implementation, testing the whole pro- 
tocol with numerical integration of the Schrodinger equa- 
tion. We show how critical observables, as the electron- 
phonon correlations can be retrieved from the trapped 
ion setup, leading to an estimation of the polaron size. 

Decomposition of the Model.- It is known that the dy- 
namics of a quantum state under the action of a Hamil- 
tonian H can be recovered by using combined fractal- 
stroboscopic symmetric decompositions [TJ [BJ. In most 
practical cases, one can assume a fractal depth of one. 
This will be the case through all the rest of our analy- 
sis. With these techniques, the target Hamiltonian H is 
decomposed in a set of m terms: H = 2<=i H m . Then, 
the symmetric decomposition for the unitary operator 
encoding the dynamics of Hamiltonian H reads 
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Here r is the degree of approximation in terms of Trotter 
steps. It has been shown [5] that, using symmetric Suzuki 
fractal decompositions, the number of gates needed to 
approximate the exact time evolution of the quantum 
state grows with the norm of the simulated Hamiltonian. 
Therefore, it is a natural problem to think of a quantum 
simulation involving particle generation, in particular of 
bosons, whose number can grow, in principle, indefinitely. 
However, in the standard approach to these problems, the 
dynamics of a bosonic Hilbert space can be recovered by 
truncating at a certain point of the number of possible 
bosonic excitations. Thus, the number of gates needed 
to achieve a certain fidelity for the simulated quantum 
state grows as more bosonic excitations are created. 
The Holstein Hamiltonian 12J, of a chain of N sites 
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Here, Ci(cj) is the annihilation (creation) operator in the 
electron site i, and bi(b\) is the phonon annihilation (cre- 
ation) operator on the site i; rii = cjcj is the electronic oc- 
cupation number operator. The parameters h, g and ujq 
stand respectively for a nearest-neighbor (NN) site hop- 
ping for the electrons, electron-phonon coupling and free 
energy of the phonons. To encode the model in a trapped- 
ion chain, we first map the fermionic operators through 
the Jordan- Wigner transformation, Cj — > YYjZi ^J ^ to 
tensor products of Pauli matrices. The mapped Hamil- 
tonian describes now a coupled spin-boson system 
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The first term can be rewritten as h 
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°f We now decompose the Hamiltonian into three 

parts, H = H± + H2 
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According to Ref. [B], one can upper bound the number 
of gates N g needed to achieve a simulation error smaller 
that e, by giving an upper bound for the norm of H |37j . 
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As mentioned before, the fractal depth k pQ can be set to 
one in most applications. Here, we show the dependence 
of the number of gates in the number of fermionic sites 
N, and on the truncation in the number of bosons M. 
As the number of created phonons increases, one needs a 
higher-level truncation, and a larger Hamiltonian norm. 
Nevertheless, this shows that we can efficiently simulate 
a 2 N x M N Hilbert space, i.e., with a number of gates 
that grows at most polynomially in N and M. To show 
the scaling of fidelities with the parameters considered, 
we plot in Fig. [T] the time dependence of the fidelity loss 
1-F(t) = l-|(* B (t)|* s (t))| 2 of the simulated wavefunc- 
tion I SI/ 5 (t)) versus the exact one \^>E(t)) as a function 
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FIG. 1. (color online), (a) Behavior of the fidelity loss 1 — 
F(t) = 1 — \{^ E(t)\^s(t)}\ 2 , for a two site configuration, as a 
function of the electron-phonon coupling strength g, for ujq = 
h/A. As the coupling increases, more phonons are created, 
the Hilbert space describing the dynamics enlarges and the 
fidelity decreases for a fixed number of approximant gates 
(r — 10 here), (b) Dependence of the fidelity loss in the 
number of sites. Here g — 0.3h, uio — 0.5h, and ten symmetric 
steps are considered (r = 10). The initial state of both plots 
corresponds to a configuration in which an electron is injected 
in the site N/2 (N even) or (N + l)/2 (N odd), and there are 
no phonons. 



of number of sites TV and coupling g. The particular de- 
composition has been chosen so that all terms in Eq. Q 
can be implemented in a linear chain of trapped ions. 

Trapped-ion setup.- We consider a set of N+l trapped 
ions in a chain, in order to simulate N fermionic sites pro- 
vided with Holstein interactions. The ions are bounded 
strongly in the radial direction, and confined longitudi- 
nally within a harmonic potential [6 . We define Ui, i — 
1, 2, ...N+l, as the frequencies of the axial normal modes. 
We relate the ion normal mode energies with the disper- 
sionless phonon energies in Eq. (|2j) via A, = — 
The three Hamiltonian steps Hi , H2 and H3 are derived 
in the interaction picture with respect to 
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where ui is the excitation energy of the individual ion 
taken as a two- level system, i.e., the carrier frequency. 
In this way, the free energies of the normal modes do not 
disappear in the interaction picture, and a flattered part 
of them is still present in order to recover the dispersion- 
less phononic spectrum. 

In order to simulate the dynamics associated to Hi 
and H 2 of Eq. Q, one has to achieve an NN Ising cou- 
pling. The possibility of obtaining an Ising field between 
linear chain of trapped ions has been proposed and re- 
alized [35]. However, in implementing NN interac- 
tions between more than two ions, one must be careful 
in designing an appropriate set of lasers and detunings 
in order to minimize the spurious non-nearest-neighbor 
(NNN) effects. To this extent, we have realized numer- 
ical simulations for a 3+1 ions setup [37], using one set 
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FIG. 2. (color online). Dynamics for the 3 + 1 ions configu- 
ration of the NN XX Hamiltonian. Dotted curves stand for 
(cj) e for the exact dynamics, and solid curves stand for {al)i 
for realistic ion interactions [i = 1,2,3 for the first, second 
and third ion). The parameters are chosen in order to have 
maxima in the fidelity F(t) = |{$ B (t)|*r(f))| 2 of ~ 0.995 
(top black curve) at time steps of ~ 333/^it. These time steps 
can be chosen as Trotter steps. 

of two pairs of counterpropagating lasers detuned close 
to the center of mass (COM) shifted mode of frequency 
Ai = v\ — ojo/3 to drive the first two ions, and another 
set of lasers detuned close to the shifted breathing mode 
of frequency A 2 , addressing the second and the third ion. 
Rabi frequencies tti of the lasers driving the i-th and the 
i + 1-th ions are chosen to achieve the desired strength 
in the Ising coupling, according to [5] , 

N-l N+l 
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In Fig. [2j the first and second ion are driven with two 
pairs of counterpropagating lasers with detuning close to 
the shifted COM mode {S 1 = ±1.0187^ for cj = h/A). 
The Rabi frequencies are chosen properly in order to 
reach a NN interaction of h/2 = O.OOl^i. Lasers driving 
the second and the third ions are detuned close to the 
shifted breathing mode at Vi — 1.731^i [6], with param- 
eters S2 = ±1.71196^1. Detunings are chosen to have a 
dynamics decoupled with respect to the phonons at time 
steps ~ 333vt and a negligible NNN interaction [37]. At 
these times, the ion spins match the exact value, phonons 
are detached from spins and the fidelity oscillation (top 
black curve) F(t) = \(^ E (t)\^ i(t))\ 2 reaches maxima, 
with peaks of ~ 0.995. 

The initial state, as in all our numerical simulations, 
is chosen to mimic a configuration in which one electron 
is injected at the center of a one dimensional lattice pro- 
vided with Holstein interactions. To this extent, all the 
spins are initialized in the opposite Z direction, except 
the one in site N/2, in case of even N, or (N + l)/2 in 
case of odd N. The spin of the last ion has to be ini- 



tialized along the Z direction in order to be a passive ion 
with respect to the dynamics, according to the protocol 
for the implementation of H3 given below. The vibra- 
tional modes are assumed to be initially cooled down to 
the ground state with resolved sideband cooling 33 . 

Notice that one can always implement a perfect NN 
coupling by using more stroboscopic steps. A possibil- 
ity is to decompose the global NN into nearest-neighbor 
pairwise interactions. Another possibility is to design a 
counter, non-nearest-neighbor interaction step between 
pairs of non-nearest neighbor ions in order to eliminate 
the spurious NNN imperfections. Given that one has an 
unwanted hijO~ % x o x , one can add more Trotter steps to 
the protocol of the form —hijo~ x a x in order to have an 
Hamiltonian free of NNN couplings. The dynamics as- 
sociated to the step with H2 is implemented similarly to 
the one of Hi , with a different choice of the initial phases 
of the lasers, in order to achieve a YY interaction. 

The Hamiltonian H3 is realized as a combination of 
2N red and blue detuned lasers with appropriate initial 
phases in order to recover a coupling of the i-th ion with 
the mj-th normal mode i r ii mi Q.ia x {b\ n .+b mi ). The i-th ion 
is driven with red and blue detuned lasers to the m^-th 
mode, establishing a one-to-one correspondence between 
the first N ions and N normal modes. Moreover, the 
last ion of the chain is driven by 2N lasers detuned in 
order to be coupled with the same modes of the ions in 
the chain. An additional rotation of the spins of all ions 
around the Y axis is applied, which can be obtained by 
acting two times with a global beam upon all the N + l 
ions at the same time. The Hamiltonian describing this 
process is, 

N 

^(^iVimA + n N+htm+hmt a^ +1 )(b mi + &tj. (8) 
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The Rabi frequencies of the lasers must be chosen accord- 
ing to tti = g/2r) itm ., fijv+M = g/^VN+i, mi - If the last 
ion is initialized with the spin aligned along the Z axis 
and not addressed by spin flip gates during the simula- 
tion, the previous described gates result in the effective 
Hamiltonian on the first N ions subspace, 

i=l 

Digital Simulation.- In general, digital protocols are 
much sensitive to the state fidelity that one can achieve 
at the end of the digital step. According to the mathe- 
matical theory, increasing the number of steps will result 
in an increased fidelity on the final simulated state. How- 
ever, if one has an error on a single step, increasing the 
number of gates will result in the accumulation of these 
errors. Thus on one hand the use of more accurate sin- 
gle gates is required, on the other hand one has to get 
a compromise between the increased fidelity due to the 
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FIG. 3. (color online). Fidelity loss for 3+1 ion configuration, 
involving Trotter simulation with perfect gates and realistic 
ion interactions, for two and three symmetric Trotter steps. 

increased number of steps and the fidelity loss due to the 
accumulated single gate error. 

To have a quantitative estimation of the fidelity loss 
with the dynamics of the full ion Hamiltonian, we have 
realized numerical integrations for the Schrodinger equa- 
tion for N = 2+1 [37] and N = 3+1 ion setups. We point 
out that we consider this reduced number of ions because 
of numerical computation restrictions, and to prove the 
feasibility of our model. In general our formalism may 
be straightforwardly extended to several ions. In Fig. 
Q, a simulation for r = 2 and r = 3 symmetric Trotter 
steps is realized. The fidelity loss 1- |(* E (i)|\E' T (t))| 2 for 
the Trotter protocol with perfect gates, i.e., associated to 
Hamiltonians Hi, H 2 and H 3 , is plotted against points 
of fidelity loss 1 - \{^ E (t)\^i(t))\ 2 obtained with realis- 
tic trapped-ion gates including the full laser interactions 
are plotted at various times. As can be appreciated, the 
fidelity loss for the ion gates is only slightly larger than 
for the exact Trotter gates, showing the feasibility of the 
protocol with realistic trapped-ion interactions. The to- 
tal simulation time has been chosen in order to remain 
under the decoherence time for the ions. The frequency of 
the center of mass mode can be assumed to be u\ ~ 2-7T x 1 
MHz. The global rotation for the ion spins can be as- 
sumed to be 7/is. The number of global rotations is 4r. 
The step for the red and blue sideband Hamiltonian can 
be performed in the same time as the step for the NN 
XX gate (or even faster). Provided with these parame- 
ters, for a final simulated time of 2000/^i ~ 318/xs, the 
time spent for the simulation can be taken of ~ 1ms. 

Tuning the coupling strength g by setting the Rabi fre- 
quencies of the red and blue detuned lasers to various val- 
ues, one can measure the different correlations between 
electron and phonon displacement at distant sites, 

X (i,j) = m)\4*(b]+b j )Mt)). (io) 

This will amounts to a signature of the polaron size [7J. 
Ranging from small to large g will lead to a measure of 



the crossover between large/small polaron. Notice that 
these correlations are mapped in our ion setup onto 

X (i,j) = (*(t)|(&^ +&J nj )^f^|*(*)), (11) 

which can be measured by mapping the motional onto the 
internal state of the auxiliary N + 1-th ion, and then de- 
tecting resonance fluorescence of ions 7V+1 and i [2"2"ll24| . 
We notice that with our setup the possibility of simulat- 
ing a 2D and 3D Holstein model is provided, by encoding 
two and three dimensional interactions into a linear chain 
by addressing distant ions with nonlocal gates @]. 

Currently, more than 100 gates have been implemented 
in a trapped-ion quantum simulation experiment with 
Trotter methods [3]. In the near future, it should be pos- 
sible to achieve hundreds or even thousands of gates per 
experiment [38 , allowing our proposal to reach about ten 
qubits. It is noteworthy to mention that our proposed 
digital quantum simulation will already overcome the 
limits of classical computers with 10 ions and 5 phonons 
per ion. In this manner, the trapped-ion quantum sim- 
ulator will prove to be a remarkable tool for simulating 
fermions coupled to bosons and related condensed-matter 
or high-energy physics scenarios. 
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SUPPLEMENTAL MATERIAL FOR "DIGITAL 
QUANTUM SIMULATION OF THE HOLSTEIN 
MODEL IN TRAPPED IONS" 

UPPER BOUND FOR THE NORM 

In this section, we give an upper bound for the the 
norm of H in Eq. (3) of the main text, in order to bound 
the error one makes with a Suzuki-Lie- Trotter expan- 
sion pQ. Consequently, we bound the number of gates 
one needs for achieving a given fidelity on the simulated 
quantum state. The norm is bounded by the sum of the 
norms of each term appearing in H. The computation of 
single norms amounts to finding the largest eigenvalue of 
the single terms 
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Let us consider the various norms separately. The 
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is bounded by the truncation in the 

number of bosons in the mode i , as is clear by the stan- 
dard Fock representation 
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(bi + b\) in the Fock basis reads 
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The characteristic polynomial of the matrix for the trun- 
cation to M bosons is given in a recursively way, 



D n (X) 
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The D n (\) are a y/2 rescaled version of the Hermite poly- 
nomials. A simple bound for the largest zero of Dm (A) 
(i.e. the norm of the bosonic displacement operator) 
is given by the expression 2y/M — 1 (see for example 
Ref. [5]). Summarizing, the norm for the Holstcin Hamil- 
tonian is upper bounded by 

\\H\\ < \h\{N- 1) + 2\g\NVM - 1 + w NM. (16) 

NUMERICAL SIMULATIONS 

In this section, we provide additional numeric plots and 
further discussions for our simulation protocol for a 2+1 
and 3+1 ion configuration. First of all, in all our numer- 
ical simulations, we have fixed the total simulated time 
to a maximum of 2000/^i, in units of the center of mass 
(COM) mode frequency v\. Assuming V\ ~ 27r X 1MHz, 
this gives a total simulated time of ~ 318^s. We choose 
Trotter steps equally extended within a time r = t/2, 
where t is the total simulated time. With these assump- 
tions, to compute the total effective simulation time, one 
has to multiply the simulated time by the number of 
terms in the decomposition of the simulated dynamics, 
i.e. 3 in our case. Then we have to add the time con- 
tribution for the global ir/4 rotations along the Y axis, 
necessary to achieve the Z-like coupling to the phonons, 
that can be estimated to be around ~ 7/xs each [3] . Con- 
sidering four global rotations per symmetric Trotter step, 
this gives a total simulation time of the order ~ 1ms for 
the r = 1 and r — 2 case. This is well below the typi- 
cal decoherence times for a trapped- ion setup [1] . Notice 
that we have made assumptions on the time extension 
for the Trotter steps, but nothing prevents to set the 
duration for the Trotter step to shorter times, as long 
as one can adjust properly the Rabi frequencies of the 
lasers used [3]. This paves the way to the scalability of 
the protocol. 

The dynamics described by Hamiltonians Hi and H 2 
can be achieved by using two pairs of counterpropagat- 
ing lasers with opposite detunings ±5j [5], driving the 
z-th and the i + 1-th ions. One can switch between a 
nearest neighbor (NN) XX/YY Ising interactions by tak- 
ing appropriate initial phases for the lasers. The effective 
spin-spin coupling generated by this kind of laser drivings 
has the form of Eq. (7) in the main text. 

In order to have negligible phonon displacements at 
the Trotter step time r, one has to choose the detuning 
5{ = ±27r/r + A m close to one the modes of (shifted) fre- 
quencies A m (thus \Si — A m | <C |A m |, \5i\). We point out 
again here that in our protocol we deal with shifted fre- 
quencies A m = v m — ujq/3, as explained in the main text, 
to take into account the desired dispersionless energies of 
the Holstein phonons. The ± sign in the choice of Si can 
be used to change the relative sign of the spin-spin inter- 
action, depending on sgn(r)i im rii+i trn ). We assume in our 
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FIG. 4. (color online). Dynamics for the 2+1 ions configu- 
ration of the NN XX Ising Hamiltonian. Dotted curves stand 
for {<jI)e for the exact dynamics, and solid curves stand for 
i for realistic ion interactions (i = 1,2 for the first and 
second ion). The parameters are chosen in order to have fi- 
delity losses of 1 - F(t) = 1 - \(9 E (t)\^i(t)}\ 2 ~ 1(T 4 (top 
black curve) at time steps of ~ 500/^i. 

simulations a relative Lamb-Dicke parameter distribution 
for ions and modes as in [6], with an overall magnitude 
of 0.1. If one chooses m — 1 for the 2+1 ions setup, i.e. a 
detuning close to the COM mode, the detuning must be 
chosen 8 = 2it/t + Ai, to obtain a positive h/2 coupling 
in the Ising NN interaction, because sgn(r]i t ir]2 } i) = + . 

In Fig. [4j we show the numerical integration for the 
dynamics of the NN XX Ising interaction for a 2+1 ions 
configuration. The simulated strength for the Ising cou- 
pling is h/2 = O.OOl^i. For uq = h/A, one has a shifted 
frequency for the COM mode of Ai = v x - 0.0005^i/3. 
By choosing r = 500/Vi, the detuning used in this case 
is S = 2twi/500 + Ai = 1.0124^, i.e. for v x = 1MHz, a 
frequency difference with the mode of A — v\ = 12.4KHz. 
The Rabi frequencies of the lasers are chosen in order to 
recover the desired strength for the Ising coupling. 

To have an idea of how the real ion interactions af- 
fect the protocol for a 2+1 ion setup, we make a plot 
of the errors on the simulated state with perfect gates 
and with ion gates in Fig. [5j One clearly sees that the 
higher fidelities obtained by using the ion gates with re- 
spect to the 3+1 ion setup are due to the higher single 
gate fidelity for the 2+1 setup, which permits to explore 
better fidelity regimes. The simulated parameters here 
are g = h/10, uj = h/4, h = 0.002^i. We remark that 
in the simulations we have used a small g/h ratio to re- 
duce the complexity of the simulation (i.e., the necessary 
truncation for the Fock space is small). Nevertheless, 
in a trapped-ion experiment, big g/h ratios with large 
freedom for the choice of ujq can be explored, thus re- 
covering the typical self trapping line for the formation 
of small polarons [7]- The time points of the ion-gate 
simulation range from t — 1000/^i to t = 2000/^i. For 
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FIG. 5. (color online). Fidelity loss for a 2+1 ions configura- 
tion as a function of time and Trotter steps. The simulated 
Holstein interaction has parameters g = h/10, ujq = h/2. Dot- 
ted and solid lines stand respectively for a simulation with one 
and two symmetric Trotter steps. Single points stand for the 
error in the simulation protocol using ion gates. 
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FIG. 6. (color online). Mean number of phonons inside a 
Trotter protocol for a 2+1 ions configuration, for the Hamil- 
tonian H = Hi + Hz , where Hi is a spin spin XX interaction 
and H$ is a spin-phonon coupling interaction, as in Eq.(4) 
in the main text. The phonons are excited within the time 
for the H3 steps (solid black line), and excited and released 
to their initial value within the Hi interactions (dotted red 
lines), with the typical oscillations for this kind of gates. 



r = 1 this gives Trotter steps ranging from r = 500/Vi to 
t = lOOO/i^i. The detuning for the NN interaction has to 
be set accordingly at each point, ranging from 1. 0124^1 
to 1.0061i/i. 

To obtain the plot of Fig. 3 in the main paper, with 
the same simulated parameter g = h/10, ujq = h/A, 
h = 0.002i/i, we have used two simultaneous NN XX 
interactions as described above. One involving the first 
two ions with a detuning close to the COM mode, and 



another one driving the second and the third ion with 
a detuning close to the breathing mode, of frequency 
V2 = 1.731i/i [6 J. To obtain an Ising interaction with 
the proper sign for the second and the third ion, we 
have to set the detuning for the second laser below 
the shifted frequency A2. This is because for a 3+1 
configuration (i.e., four ions in a linear trap), one has 
sgn[rjuri2i\ ^ sgn[r]2, 2^32] |6]. The detuning for the sec- 
ond set of lasers is therefore chosen to be 5 = — 27r/r+A 2 . 
For example, for the simulation point at t = 1000/t/i, cor- 
responding to t — 250/^i for r — 2, one has 5% — 1.025j/i 
and 82 = 1.7057^i. Using these parameters, it turns out 
that the spurious non-nearest-neighbor coupling between 
the first and the third ion is negligible. 



To get an insight of what happens to the phonon 
population of the COM mode inside a Trotter protocol, 
one can have a look to Fig. [6| Here, it is shown the 
mean number of phonons for a 2+1 ion setup using a 
symmetric decomposition at r = 1 of the Hamiltonian 
H = Hi + H3, where Hi is a XX Ising interaction ob- 
tained with a detuning close to Ai and H% is a Z-like cou- 
pling to phonons. The decomposed evolution operator 
has the form U 2 (t) = e -iHit/2 e -iH 3 t/2 e -iH 3 t/2 e -iHtt/2^ 

We see that, in the first and the last step, the two Ising 
interactions create phonons, while relaxing them at the 
end of the step, because the laser detuning and Rabi 
frequencies are chosen to obtain detachment from the 
phonons at the end of the Trotter step. In the two mid- 
dle steps, the phonons are excited according to the #3 
Hamiltonian. The final mean value for the phonon num- 
ber is recovered with respect to the dashed line value, 
which is the numerical value according to the exact evo- 
lution operator e~ lHt , with an error of ~ 1%. Notice 
that since the decomposition involves symmetric Trotter 
steps, and each one is chosen to be of the same duration, 
the total simulation time is doubled. 
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